import pandas as pd
import geopandas as gpd
import rasterio
import numpy as np
from shapely.geometry import Point
import libpysal as lp
from mgwr.gwr import GWR
import matplotlib.pyplot as plt

# 读取表格数据
df = pd.read_excel(r'D:\data1\processed_data.xlsx')

# 提取坐标和蓄积量
geometry = [Point(xy) for xy in zip(df['纵坐标'], df['横坐标'])]
gdf = gpd.GeoDataFrame(df, geometry=geometry)

# 读取NDVI栅格数据
ndvi_raster = rasterio.open(r'D:\data1\hunan_ndvi_2021.tif')

# 提取每个样地的NDVI值
ndvi_values = []
for idx, row in gdf.iterrows():
    lon, lat = row['纵坐标'], row['横坐标']
    try:
        # 转换为栅格坐标
        row_idx, col_idx = ndvi_raster.index(lon, lat)
        ndvi = ndvi_raster.read(1)[row_idx, col_idx]
        ndvi_values.append(ndvi)
    except:
        ndvi_values.append(np.nan)

# 将NDVI值添加到GeoDataFrame
gdf['NDVI'] = ndvi_values

# 筛选有效数据（去除NaN值）
gdf = gdf.dropna(subset=['NDVI', '蓄积量'])

# 提取因变量和自变量
y = gdf['蓄积量'].values.reshape(-1, 1)
X = gdf[['NDVI']].values

# 计算坐标矩阵
coords = np.array(list(zip(gdf.geometry.x, gdf.geometry.y)))

# 构建空间权重矩阵
bw = 50  # 初始带宽，可以使用bw = gwr.sel_bw来自动选择
gwr_model = GWR(coords, y, X, bw=bw, family='gaussian', kernel='bisquare')
gwr_results = gwr_model.fit()

# 输出模型结果
print("GWR模型结果:")
print(f"R²: {gwr_results.R2}")
print(f"AIC: {gwr_results.AIC}")

# 预测值
gdf['Predicted'] = gwr_results.predictions.flatten()

# 可视化预测结果
plt.figure(figsize=(12, 6))

# 散点图：预测值 vs 实际值
plt.subplot(1, 2, 1)
plt.scatter(gdf['蓄积量'], gdf['Predicted'], alpha=0.5)
plt.xlabel('实际蓄积量')
plt.ylabel('预测蓄积量')
plt.title('预测值 vs 实际值')

# 残差图
plt.subplot(1, 2, 2)
residuals = gdf['蓄积量'] - gdf['Predicted']
plt.scatter(gdf.geometry.x, residuals, alpha=0.5)
plt.xlabel('经度')
plt.ylabel('残差')
plt.title('残差空间分布')

plt.tight_layout()
plt.show()